1 Load libraries

# load libraries
library(tictoc)
library(tidyverse)
library(terra)
library(sf)
library(RColorBrewer)
library(rts)
library(doParallel)
library(parallel)
library(foreach)
library(data.table)

2 Functions

# Function to save layers
saveLayers = function(dataRast, shapeFile, pathFile, strRemove,
                      x_min=xmin(dataRast), x_max=xmax(dataRast), 
                      y_min=ymin(dataRast), y_max=ymax(dataRast)){
  
  # Iterate all the layers
  for (i in 1:nlyr(dataRast)){
    print(i)
    dataRast.i <- dataRast[[i]]
    # Crop data
    dataRast.basin.i <- mask(dataRast.i, mask = shapeFile) %>%
      crop(ext(x_min, x_max, y_min, y_max))
    # Path of the file
    path <- gsub(strRemove, "", names(dataRast.basin.i)) %>%
      paste0(pathFile, ., ".tif")
    # Save file
    writeRaster(dataRast.basin.i, path, overwrite=TRUE)
    # Clear memory
    gc()
  }
}

# Function to rename layers
renameLayers <- function(dataRast, fileStart, prefix){
  # index of file name in the path
  id.date <- unlist(gregexpr(fileStart, sources(dataRast)[[1]])) + nchar(fileStart)
  # rename layers
  names(dataRast) <- 
    substr(sources(dataRast), id.date, (id.date+50)) %>%
    gsub(".tif", "_1", .) %>%
    as.Date("%Y_%m_%d") %>%
    format(., '%Y_%m') %>%
    paste0(prefix, .)
  
  return(dataRast)
}

3 Load data

load("R_data/data-amaz_na3.RData")

#=====================

4 Import shape file

# Import shape file
amaz.basin.shp <- st_read("Amazon_new_data/0. Amazon_shapefile/projected/amazon_shp_projected.shp")
## Reading layer `amazon_shp_projected' from data source 
##   `/Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/0. Amazon_shapefile/projected/amazon_shp_projected.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -2156811 ymin: 1625314 xmax: 1745999 ymax: 4555427
## Projected CRS: South_America_Albers_Equal_Area_Conic

#=====================

5 Burnt Area data

5.1 Import data

# list of files
amaz.burntArea.list <- list.files("Amazon_new_data/1. Burnt Area/03. Working Data",
                                  full.names=TRUE,
                                  pattern = ".tif$")
# Import data with "Terra"
burntArea.rast <- rast(amaz.burntArea.list)
burntArea.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : burntarea_working_2001_1.tif  
##               burntarea_working_2001_10.tif  
##               burntarea_working_2001_11.tif  
##               ... and 237 more source(s)
## names       : fire_~_proj, fire_~_proj, fire_~_proj, fire_~_proj, fire_~_proj, fire_~_proj, ... 
## min values  :          -2,          -2,          -2,          -2,          -2,          -2, ... 
## max values  :           1,           1,           1,           1,           1,           1, ...

5.2 Verification of the data

5.2.1 Names of layers

names(burntArea.rast) %>% head()
## [1] "fire_2001_05_06_1_proj"  "fire_2001_05_06_10_proj"
## [3] "fire_2001_05_06_11_proj" "fire_2001_05_06_12_proj"
## [5] "fire_2001_05_06_2_proj"  "fire_2001_05_06_3_proj"

5.2.2 Find duplicate name layers

idx.dplc <- names(burntArea.rast) %>% duplicated() %>% which()
idx.dplc
## [1] 144

5.2.4 Verification of the values

# Verification of the values
burntArea.minmax <- minmax(burntArea.rast) %>% t() %>% as.data.frame()
burntArea.minmax
burntArea.minmax[which((burntArea.minmax[,1] != -2) & (burntArea.minmax[,2] != 1)),]

5.2.5 Plot of the months September and November 2012

# Create palette of color
pal <- colorRampPalette(c("lightblue", "forestgreen", "firebrick"))
# burntArea.freq <- freq(burntArea.rast, digits=0, usenames=T)
burntArea.freq[(burntArea.freq[,2] != -2)&(burntArea.freq[,2] != 0)&(burntArea.freq[,2] != 1),]
plot(burntArea.rast[[c(142,144)]], col=pal(3))

5.3 Modification of the data: Rename layers

5.3.1 Rename layers

# Rename layers
burntArea.rast <- renameLayers(burntArea.rast, 'burntarea_working_', 'fire_')
# Order layers
ordered.names.1 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>%
  paste0("fire_", .)
burntArea.rast <- burntArea.rast[[ordered.names.1]]
burntArea.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : burntarea_working_2001_1.tif  
##               burntarea_working_2001_2.tif  
##               burntarea_working_2001_3.tif  
##               ... and 237 more source(s)
## names       : fire_2001_01, fire_2001_02, fire_2001_03, fire_2001_04, fire_2001_05, fire_2001_06, ... 
## min values  :           -2,           -2,           -2,           -2,           -2,           -2, ... 
## max values  :            1,            1,            1,            1,            1,            1, ...

5.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
burntArea.rts <- rts(burntArea.rast, seq.dates)
burntArea.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/1. Burnt Area/03. Working Data/burntarea_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/1. Burnt Area/03. Working Data/burntarea_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/1. Burnt Area/03. Working Data/burntarea_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : -2 -2 -2 -2 -2 -2 -2 -2 -2 -2 ... 
## max values  :  1  1  1  1  1  1  1  1  1  1 ...

5.4 Plot

5.4.1 Define palette of color

# Create palette of color
pal <- colorRampPalette(c("lightblue", "forestgreen", "firebrick"))
pal2 <- colorRampPalette(c("forestgreen", "firebrick"))

5.4.2 Plot all the area

# Plot Amazon in the year 2012.
plot(burntArea.rts[['2012']], col=pal(3))

# plot Amazon in October 2020.
plot(burntArea.rts[['2020-10-01']], col=pal(3), main="Burnt area in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
burntArea.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=pal(3), main="Burnt area in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

5.4.3 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
burntArea.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=pal(3))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
burntArea.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=pal(3), main="Burnt area in October 2020")
plot(amaz.basin.shp$geometry, add=T)

5.5 Repalce undifined values by NA

burntArea.rast[[c(139, 141)]][burntArea.rast[[c(139, 141)]] == -1] <- NA
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          
burntArea.rast[[c(139, 141)]]
## class       : SpatRaster 
## dimensions  : 5860, 7806, 2  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : spat_r1U8IiItX34kjKM_10135.tif  
##               spat_r1U8IiItX34kjKM_10135.tif  
## names       : fire_2012_07, fire_2012_09 
## min values  :          NaN,          NaN 
## max values  :          NaN,          NaN

5.6 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("burntArea")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.burntArea.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'burntarea_working_', '')
#   # Replace negative value by `NA`
#   if (names(ras) %in% c("2012_07", "2012_09")) {ras[ras == -1] <- NA}
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# burntArea.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(burntArea.freq.na)[3] <- "burntArea_na"
# burntArea.freq.na <- burntArea.freq.na[order(burntArea.freq.na$layer)]
burntArea.freq.na

#=====================

6 Land Cover data

6.1 Import data

# list of files
amaz.landCover.list <- list.files("Amazon_new_data/2. Land Cover/03. Working Data",
                                  full.names=TRUE,
                                  pattern = ".tif$")
# Import data with "Terra"
landCover.rast <- rast(amaz.landCover.list)
landCover.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : landcover_working_2001_1.tif  
##               landcover_working_2001_10.tif  
##               landcover_working_2001_11.tif  
##               ... and 237 more source(s)
## names       : lulc_~_proj, lulc_~_proj, lulc_~_proj, lulc_~_proj, lulc_~_proj, lulc_~_proj, ... 
## min values  :           0,           0,           0,           0,           0,           0, ... 
## max values  :          10,          10,          10,          10,          10,          10, ...

6.2 Verification of the data

6.2.1 Verification the Names of layers

names(landCover.rast) %>% head(24)
##  [1] "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj"
##  [5] "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj"
##  [9] "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj" "lulc_2001_proj"
## [13] "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj"
## [17] "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj"
## [21] "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj" "lulc_2002_proj"

6.2.2 Verification of the values

# Verification of the values
landCover.minmax <- minmax(landCover.rast) %>% t() %>% as.data.frame()
landCover.minmax

6.2.3 Frequency of each value

# landCover.freq <- freq(landCover.rast, digits=0, usenames=T)
landCover.freq
# 
landCover.freq[which(!(landCover.freq$value %in% c(0,1,2,3,4,5,6,7,8,9,10)) ),]

6.3 Modification of the data

6.3.1 Rename layers

# Rename layers
landCover.rast <- renameLayers(landCover.rast, 'landcover_working_', 'lulc_')
# Order layers
ordered.names.2 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("lulc_", .)
landCover.rast <- landCover.rast[[ordered.names.2]]
landCover.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : landcover_working_2001_1.tif  
##               landcover_working_2001_2.tif  
##               landcover_working_2001_3.tif  
##               ... and 237 more source(s)
## names       : lulc_2001_01, lulc_2001_02, lulc_2001_03, lulc_2001_04, lulc_2001_05, lulc_2001_06, ... 
## min values  :            0,            0,            0,            0,            0,            0, ... 
## max values  :           10,           10,           10,           10,           10,           10, ...

6.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
landCover.rts <- rts(landCover.rast, seq.dates)
landCover.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/2. Land Cover/03. Working Data/landcover_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/2. Land Cover/03. Working Data/landcover_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/2. Land Cover/03. Working Data/landcover_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : 0 0 0 0 0 0 0 0 0 0 ... 
## max values  : 10 10 10 10 10 10 10 10 10 10 ...

6.4 Plot

6.4.1 Plot all the area

# Plot Amazon in the year 2020
plot(landCover.rts[['2020']], col=colorRamps::blue2green(10))

# plot Amazon in October 2020.
plot(landCover.rts[['2020-10-01']], col=colorRamps::blue2green(10), main="Land cover in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
landCover.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=colorRamps::blue2green(10), main="Land cover in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

6.4.2 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
landCover.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::blue2green(10))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
landCover.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::blue2green(10), main="Land cover in October 2020")
plot(amaz.basin.shp$geometry, add=T)

6.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("landCover")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.landCover.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'landcover_working_', '')
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# landCover.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(landCover.freq.na)[3] <- "landCover_na"
# landCover.freq.na <- landCover.freq.na[order(landCover.freq.na$layer)]
landCover.freq.na

#=====================

7 Precipitation data

7.1 Import data

# list of files
amaz.precipitation.list <- list.files("Amazon_new_data/3. Precipitation/03. Working Data",
                                      full.names=TRUE,
                                      pattern = ".tif$")
# Import data with "Terra"
precipitation.rast <- rast(amaz.precipitation.list)
precipitation.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : precipitation_working_2001_1.tif  
##               precipitation_working_2001_10.tif  
##               precipitation_working_2001_11.tif  
##               ... and 237 more source(s)
## names       : prec_~_proj, prec_~_proj, prec_~_proj, prec_~_proj, prec_~_proj, prec_~_proj, ... 
## min values  :           0,           0,           0,           0,           0,           0, ... 
## max values  :        1461,        1735,        1828,        1934,        1433,        1390, ...

7.2 Verification of the data

7.2.1 Verification the Names of layers

names(precipitation.rast) %>% head()
## [1] "prec_2001_1_proj"  "prec_2001_10_proj" "prec_2001_11_proj"
## [4] "prec_2001_12_proj" "prec_2001_2_proj"  "prec_2001_3_proj"

7.2.2 Verification of the values

# Verification of the values
precipitation.minmax <- minmax(precipitation.rast) %>% t() %>% as.data.frame()
precipitation.minmax

7.3 Modification of the data

7.3.1 Rename layers

# Rename layers
precipitation.rast <- renameLayers(precipitation.rast, 'precipitation_working_', 'prec_')
# Order layers
ordered.names.3 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("prec_", .)
precipitation.rast <- precipitation.rast[[ordered.names.3]]
precipitation.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : precipitation_working_2001_1.tif  
##               precipitation_working_2001_2.tif  
##               precipitation_working_2001_3.tif  
##               ... and 237 more source(s)
## names       : prec_2001_01, prec_2001_02, prec_2001_03, prec_2001_04, prec_2001_05, prec_2001_06, ... 
## min values  :            0,            0,            0,            0,            0,            0, ... 
## max values  :         1461,         1433,         1390,         1096,         1678,         1706, ...

7.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
precipitation.rts <- rts(precipitation.rast, seq.dates)
precipitation.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : 0 0 0 0 0 0 0 0 0 0 ... 
## max values  : 1461 1433 1390 1096 1678 1706 1807 1988 2230 1735 ...

7.4 Plot

7.4.1 Plot all the area

# Plot Amazon in the year 2020
plot(precipitation.rts[['2020']], col=grDevices::cm.colors(1000))

# plot Amazon in October 2020.
plot(precipitation.rts[['2020-10-01']], col=grDevices::cm.colors(1000), main="Precipitation in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
precipitation.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=grDevices::cm.colors(1000), main="Precipitation in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

7.4.2 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
precipitation.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=grDevices::cm.colors(1000))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
precipitation.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=grDevices::cm.colors(1000), main="Precipitation in October 2020")
plot(amaz.basin.shp$geometry, add=T)

7.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("precipitation")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.precipitation.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'precipitation_working_', '')
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# precipitation.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(precipitation.freq.na)[3] <- "precipitation_na"
# precipitation.freq.na <- precipitation.freq.na[order(precipitation.freq.na$layer)]
precipitation.freq.na

7.6 Sample of NA values (February 2002)

7.6.1 Compute

precipitation.2002.02.name <- c("Amazon_new_data/3. Precipitation/03. Working Data/precipitation_working_2002_2.tif")
precipitation.2002.02.rast <- rast(precipitation.2002.02.name)
precipitation.2002.02.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : precipitation_working_2002_2.tif 
## name        : prec_2002_2_proj 
## min value   :                0 
## max value   :             1439
#
precipitation.2002.02.nonNA <- not.na(precipitation.2002.02.rast)
## 
|---------|---------|---------|---------|
=========================================
                                          
precipitation.2002.02.nonNA.mask <- mask(precipitation.2002.02.nonNA, amaz.basin.shp)
## 
|---------|---------|---------|---------|
=========================================
                                          
#
precipitation.2002.02.nonNA.df <- as.data.frame(precipitation.2002.02.nonNA.mask, xy=T)
precipitation.2002.02.NA.df <- precipitation.2002.02.nonNA.df[precipitation.2002.02.nonNA.df[,3]==F,]
precipitation.2002.02.NA.df
dim(precipitation.2002.02.NA.df)[1]
## [1] 34

7.6.2 Plot

# Plot
plot(precipitation.2002.02.rast)

plot(precipitation.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
plot(amaz.basin.shp$geometry, add=T)

#
x_min=-0.45e+06; x_max=-0.3e+06; y_min=ymin(precipitation.2002.02.nonNA.mask); y_max=1.71e+06
precipitation.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#=====================

8 Soil Moisture data

8.1 Import data

# list of files
amaz.soilMoisture.list <- list.files("Amazon_new_data/4. Soil Moisture/03. Working Data",
                                     full.names=TRUE,
                                     pattern = ".tif$")
# Import data with "Terra"
soilMoisture.rast <- rast(amaz.soilMoisture.list)
soilMoisture.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : soilmoisture_working_2001_1.tif  
##               soilmoisture_working_2001_10.tif  
##               soilmoisture_working_2001_11.tif  
##               ... and 237 more source(s)
## names       :   soilm~_proj,   soilm~_proj,   soilm~_proj,   soilm~_proj,   soilm~_proj,   soilm~_proj, ... 
## min values  : -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, ... 
## max values  :  1.076058e+03,  4.290171e+03,  7.905507e+02,  5.636223e+02,  7.509722e+02,  6.363759e+02, ...

8.2 Verification of the data

8.2.1 Verification the Names of layers

names(soilMoisture.rast) %>% head()
## [1] "soilm_2001_1_proj"  "soilm_2001_10_proj" "soilm_2001_11_proj"
## [4] "soilm_2001_12_proj" "soilm_2001_2_proj"  "soilm_2001_3_proj"

8.2.2 Verification of the values

# Verification of the values
soilMoisture.minmax <- minmax(soilMoisture.rast) %>% t() %>% as.data.frame()
soilMoisture.minmax

8.2.3 Frequency

# soilMoisture.freq <- freq(soilMoisture.rast, digits=3, usenames=T)
soilMoisture.freq
#
soilMoisture.freq[soilMoisture.freq$value < 0,]

8.3 Modification of the data

8.3.1 Rename layers

# Rename layers
soilMoisture.rast <- renameLayers(soilMoisture.rast, 'soilmoisture_working_', 'soilm_')
# Order layers
ordered.names.4 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("soilm_", .)
soilMoisture.rast <- soilMoisture.rast[[ordered.names.4]]
# Print
soilMoisture.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : soilmoisture_working_2001_1.tif  
##               soilmoisture_working_2001_2.tif  
##               soilmoisture_working_2001_3.tif  
##               ... and 237 more source(s)
## names       :   soilm~01_01,   soilm~01_02,   soilm~01_03,   soilm~01_04,   soilm~01_05,   soilm~01_06, ... 
## min values  : -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, -9.990000e+08, ... 
## max values  :  1.076058e+03,  7.509722e+02,  6.363759e+02,  7.215025e+02,  6.892701e+02,  5.409775e+02, ...
cat('\n---------- \n\n')
## 
## ----------
names(soilMoisture.rast) %>% head()
## [1] "soilm_2001_01" "soilm_2001_02" "soilm_2001_03" "soilm_2001_04"
## [5] "soilm_2001_05" "soilm_2001_06"

8.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
soilMoisture.rts <- rts(soilMoisture.rast, seq.dates)
soilMoisture.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 -1e+09 ... 
## max values  : 1076  751  636  722  689  541  713 1167 1568 4290 ...

8.4 Plot

8.4.1 Plot all the area

# Plot Amazon in the year 2020
plot(soilMoisture.rts[['2020']], col=colorRamps::matlab.like2(1000))

# plot Amazon in October 2020.
plot(soilMoisture.rts[['2020-10-01']], col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
soilMoisture.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

#
x_min=0.1e+06; x_max=0.5e+06; y_min=4.25e+06; y_max=ymax(burntArea.rast)
soilMoisture.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

8.4.2 Repalce negative values by NA for 2020

soilMoisture.2020.rts <- soilMoisture.rts[['2020']]
# Remove negative values.
soilMoisture.2020.rts@raster[soilMoisture.2020.rts@raster < 0] <- NA
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          
soilMoisture.2020.rts
## Raster Time Series with monthly periodicity from 2020-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /private/var/folders/91/ckb_ktpx5mzcsls1l4tmm5xm0000gp/T/Rtmp6sflAg/spat_CKF46jJIMXAcKYP_10135.tif 
## raster dimensions  : 5860, 7806, 12  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : 0.00645 0.00486 0.00368 0.00268 0.00203 0.00159 0.00128 0.00104 0.00084 0.00067 ... 
## max values  : 236 379 443 413 777 431 482 494 544 258 ...

8.4.3 Plot all the area

# Plot Amazon in the year 2020
plot(soilMoisture.2020.rts, col=colorRamps::matlab.like2(1000))

# plot Amazon in October 2020.
plot(soilMoisture.2020.rts[['2020-10-01']], col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
soilMoisture.2020.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

8.4.4 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
soilMoisture.2020.rts@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::matlab.like2(1000))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
soilMoisture.2020.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::matlab.like2(1000), main="Soil Moisture in October 2020")
plot(amaz.basin.shp$geometry, add=T)

8.5 Repalce negative values by NA for all data

# # Remove negative values.
# soilMoisture.rast[soilMoisture.rast < 0] <- NA
# # Save raster file
# writeRaster(soilMoisture.rast, "Amazon_new_data/4. Soil Moisture/soilmoisture_working_na.tif", overwrite=TRUE)

# Load data
soilMoisture.na.rast <- rast("Amazon_new_data/4. Soil Moisture/soilmoisture_working_na.tif")
soilMoisture.na.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : soilmoisture_working_na.tif 
## names       :  soilm~01_01,  soilm~01_02,  soilm~01_03,  soilm~01_04, soilm~01_05, soilm~01_06, ... 
## min values  : 3.900412e-02,   0.05741884,   0.07022829,   0.05813476,   0.1009746,   0.1354199, ... 
## max values  : 1.076058e+03, 750.97216797, 636.37585449, 721.50250244, 689.2701416, 540.9775391, ...

8.6 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("soilmoisture")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.soilMoisture.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'soilmoisture_working_', '')
#   # Replace negative value by `NA`
#   ras[ras < 0] <- NA
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# soilmoisture.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(soilmoisture.freq.na)[3] <- "soilmoisture_na"
# soilmoisture.freq.na <- soilmoisture.freq.na[order(soilmoisture.freq.na$layer)]
soilmoisture.freq.na

8.7 Sample of NA values (February 2002)

8.7.1 Compute

soilmoisture.2002.02.name <- c("Amazon_new_data/4. Soil Moisture/03. Working Data/soilmoisture_working_2002_2.tif")
soilmoisture.2002.02.rast <- rast(soilmoisture.2002.02.name)
# Repalce negative values by `NA` for all data
soilmoisture.2002.02.rast[soilmoisture.2002.02.rast < 0] <- NA
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          
soilmoisture.2002.02.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : spat_8KHrIMKbEZfP8DW_10135.tif 
## name        : soilm_2002_2_proj 
## min value   :        0.06081576 
## max value   :      465.46282959
#
soilmoisture.2002.02.nonNA <- not.na(soilmoisture.2002.02.rast)
## 
|---------|---------|---------|---------|
=========================================
                                          
soilmoisture.2002.02.nonNA.mask <- mask(soilmoisture.2002.02.nonNA, amaz.basin.shp)
## 
|---------|---------|---------|---------|
=========================================
                                          
#
soilmoisture.2002.02.nonNA.df <- as.data.frame(soilmoisture.2002.02.nonNA.mask, xy=T)
soilmoisture.2002.02.NA.df <- soilmoisture.2002.02.nonNA.df[soilmoisture.2002.02.nonNA.df[,3]==F,]
soilmoisture.2002.02.NA.df
dim(soilmoisture.2002.02.NA.df)[1]
## [1] 1975

8.7.2 Plot

# Plot
plot(soilmoisture.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))

#
x_min=-0.45e+06; x_max=-0.3e+06; y_min=ymin(soilmoisture.2002.02.nonNA.mask); y_max=1.71e+06
soilmoisture.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#
x_min=1.2e+06; x_max=2.5e+06; y_min=3e+06; y_max=3.5e+06
soilmoisture.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#
x_min=0.1e+06; x_max=0.5e+06; y_min=4.25e+06; y_max=ymax(soilmoisture.2002.02.nonNA.mask)
soilmoisture.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#=====================

9 Elevation data

9.1 Import data

# list of files
amaz.elevation.list <- list.files("Amazon_new_data/5. Elevation/03. Working Data",
                                  full.names=TRUE,
                                  pattern = ".tif$")
# Import data with "Terra"
elevation.rast <- rast(amaz.elevation.list)
elevation.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : elevation_working.tif 
## name        : elevation_proj 
## min value   :          -85.5 
## max value   :         6470.5

9.2 Verification of the data

9.2.1 Verification the Names of layers

names(elevation.rast)
## [1] "elevation_proj"

9.2.2 Verification of the values

# Verification of the values
elevation.minmax <- minmax(elevation.rast) %>% t() %>% as.data.frame()
elevation.minmax

9.3 Modification of the data

9.3.1 Rename layers

names(elevation.rast) <- "elevation"
elevation.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : elevation_working.tif 
## name        : elevation 
## min value   :     -85.5 
## max value   :    6470.5

9.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2020-12-1"), as.Date("2020-12-1"), by = "month")
elevation.rts <- rts(elevation.rast, seq.dates)
elevation.rts
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/5. Elevation/03. Working Data/elevation_working.tif 
## raster dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : -86 
## max values  : 6470

9.4 Plot

9.4.1 Plot all the area

# plot.
plot(elevation.rts[['2020-12-01']], col=colorRamps::matlab.like(1000), main="Elevation")
plot(amaz.basin.shp$geometry, add=T)

#
elevation.rts[['2020-12-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=colorRamps::matlab.like(1000), main="Elevation")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

9.4.2 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon.
elevation.rts[['2020-12-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::matlab.like(1000), main="Elevation")
plot(amaz.basin.shp$geometry, add=T)

9.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("elevation")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.elevation.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   names(ras) <- "2020_12"
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# elevation.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(elevation.freq.na)[3] <- "elevation_na"
# elevation.freq.na <- elevation.freq.na[order(elevation.freq.na$layer)]
elevation.freq.na

#=====================

10 LandSurfaceTemp data

10.1 Import data

# list of files
amaz.landSurfaceTemp.list <- list.files("Amazon_new_data/6. LandSurfaceTemp/03. Working Data",
                                     full.names=TRUE,
                                     pattern = ".tif$")
# Import data with "Terra"
landSurfaceTemp.rast <- rast(amaz.landSurfaceTemp.list)
landSurfaceTemp.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : landsurftemp_working_2001_1.tif  
##               landsurftemp_working_2001_10.tif  
##               landsurftemp_working_2001_11.tif  
##               ... and 237 more source(s)
## names       : Month~ature, Month~ature, Month~ature, Month~ature, Month~ature, Month~ature, ... 
## min values  :       13255,       13746,       13376,       13395,       12869,       13508, ... 
## max values  :       16387,       16415,       16356,       16434,       16432,       16455, ...

10.2 Verification of the data

10.2.1 Verification the Names of layers

names(landSurfaceTemp.rast) %>% head()
## [1] "Monthly daytime 3min CMG Land-surface Temperature"
## [2] "Monthly daytime 3min CMG Land-surface Temperature"
## [3] "Monthly daytime 3min CMG Land-surface Temperature"
## [4] "Monthly daytime 3min CMG Land-surface Temperature"
## [5] "Monthly daytime 3min CMG Land-surface Temperature"
## [6] "Monthly daytime 3min CMG Land-surface Temperature"

10.2.2 Verification of the values

# Verification of the values
landSurfaceTemp.minmax <- minmax(landSurfaceTemp.rast) %>% t() %>% as.data.frame()
landSurfaceTemp.minmax

10.3 Modification of the data

10.3.1 Rename layers

landSurfaceTemp.rast <- renameLayers(landSurfaceTemp.rast, 'landsurftemp_working_', 'landsurftemp_')
# Order layers
ordered.names.6 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("landsurftemp_", .)
landSurfaceTemp.rast <- landSurfaceTemp.rast[[ordered.names.6]]
# Print
landSurfaceTemp.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : landsurftemp_working_2001_1.tif  
##               landsurftemp_working_2001_2.tif  
##               landsurftemp_working_2001_3.tif  
##               ... and 237 more source(s)
## names       : lands~01_01, lands~01_02, lands~01_03, lands~01_04, lands~01_05, lands~01_06, ... 
## min values  :       13255,       12869,       13508,       13464,       13727,       13364, ... 
## max values  :       16387,       16432,       16455,       16223,       15990,       15897, ...
cat('\n---------- \n\n')
## 
## ----------
names(landSurfaceTemp.rast) %>% head()
## [1] "landsurftemp_2001_01" "landsurftemp_2001_02" "landsurftemp_2001_03"
## [4] "landsurftemp_2001_04" "landsurftemp_2001_05" "landsurftemp_2001_06"

10.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
landSurfaceTemp.rts <- rts(landSurfaceTemp.rast, seq.dates)
landSurfaceTemp.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : 13255 12869 13508 13464 13727 13364 13593 13570 13800 13746 ... 
## max values  : 16387 16432 16455 16223 15990 15897 15908 16023 16139 16415 ...

10.4 Plot

10.4.1 Plot all the area

# Plot Amazon in the year 2020
plot(landSurfaceTemp.rts[['2020']], col=colorRamps::green2red(1000))

# plot Amazon in October 2020.
plot(landSurfaceTemp.rts[['2020-10-01']], col=colorRamps::green2red(1000), main="Land Surface Temp in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
landSurfaceTemp.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=colorRamps::green2red(1000), main="Land Surface Temp in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

10.4.2 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
landSurfaceTemp.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::green2red(1000))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
landSurfaceTemp.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::green2red(1000), main="Land Surface Temp in October 2020")
plot(amaz.basin.shp$geometry, add=T)

10.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("landsurftemp")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.landSurfaceTemp.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'landsurftemp_working_', '')
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# landsurftemp.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(landsurftemp.freq.na)[3] <- "landsurftemp_na"
# landsurftemp.freq.na <- landsurftemp.freq.na[order(landsurftemp.freq.na$layer)]
landsurftemp.freq.na

10.6 Sample of NA values (February 2002)

10.6.1 Compute

landSurfaceTemp.2002.02.name <- c("Amazon_new_data/6. LandSurfaceTemp/03. Working Data/landsurftemp_working_2002_2.tif")
landSurfaceTemp.2002.02.rast <- rast(landSurfaceTemp.2002.02.name)
landSurfaceTemp.2002.02.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : landsurftemp_working_2002_2.tif 
## name        : Monthly daytime 3min CMG Land-surface Temperature 
## min value   :                                             13421 
## max value   :                                             16484
#
landSurfaceTemp.2002.02.nonNA <- not.na(landSurfaceTemp.2002.02.rast)
## 
|---------|---------|---------|---------|
=========================================
                                          
landSurfaceTemp.2002.02.nonNA.mask <- mask(landSurfaceTemp.2002.02.nonNA, amaz.basin.shp)
## 
|---------|---------|---------|---------|
=========================================
                                          
#
landSurfaceTemp.2002.02.nonNA.df <- as.data.frame(landSurfaceTemp.2002.02.nonNA.mask, xy=T)
landSurfaceTemp.2002.02.NA.df <- landSurfaceTemp.2002.02.nonNA.df[landSurfaceTemp.2002.02.nonNA.df[,3]==F,]
landSurfaceTemp.2002.02.NA.df
dim(landSurfaceTemp.2002.02.NA.df)[1]
## [1] 2085188

10.6.2 Plot

# Plot
plot(landSurfaceTemp.2002.02.rast)

plot(landSurfaceTemp.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))
plot(amaz.basin.shp$geometry, add=T)

#=====================

11 Specific Humidity data

11.1 Import data

# list of files
amaz.humidity.list <- list.files("Amazon_new_data/7. Specific Humidity/03. Working Data",
                                 full.names=TRUE,
                                 pattern = ".tif$")
# Import data with "Terra"
humidity.rast <- rast(amaz.humidity.list)
humidity.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : humidity_working_2001_1.tif  
##               humidity_working_2001_10.tif  
##               humidity_working_2001_11.tif  
##               ... and 237 more source(s)
## names       : humid~_proj, humid~_proj, humid~_proj, humid~_proj, humid~_proj, humid~_proj, ... 
## min values  : 0.003805317, 0.002230128, 0.002613974,  0.00270508, 0.004156957, 0.004153194, ... 
## max values  : 0.018998224, 0.019306520, 0.019978305,  0.02033914, 0.018878255, 0.019174447, ...

11.2 Verification of the data

11.2.1 Verification the Names of layers

names(humidity.rast) %>% head()
## [1] "humidity_2001_1_proj"  "humidity_2001_10_proj" "humidity_2001_11_proj"
## [4] "humidity_2001_12_proj" "humidity_2001_2_proj"  "humidity_2001_3_proj"

11.2.2 Verification of the values

# Verification of the values
humidity.minmax <- minmax(humidity.rast) %>% t() %>% as.data.frame()
humidity.minmax

11.3 Modification of the data

11.3.1 Rename layers

humidity.rast <- renameLayers(humidity.rast, 'humidity_working_', 'humidity_')
# Order layers
ordered.names.7 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("humidity_", .)
humidity.rast <- humidity.rast[[ordered.names.7]]
# Print
humidity.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : humidity_working_2001_1.tif  
##               humidity_working_2001_2.tif  
##               humidity_working_2001_3.tif  
##               ... and 237 more source(s)
## names       : humid~01_01, humid~01_02, humid~01_03, humid~01_04, humid~01_05, humid~01_06, ... 
## min values  : 0.003805317, 0.004156957, 0.004153194, 0.003362617, 0.002170095, 0.001819314, ... 
## max values  : 0.018998224, 0.018878255, 0.019174447, 0.020020029, 0.019958658, 0.019542987, ...
cat('\n---------- \n\n')
## 
## ----------
names(humidity.rast) %>% head()
## [1] "humidity_2001_01" "humidity_2001_02" "humidity_2001_03" "humidity_2001_04"
## [5] "humidity_2001_05" "humidity_2001_06"

11.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
humidity.rts <- rts(humidity.rast, seq.dates)
humidity.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : 0.00381 0.00416 0.00415 0.00336 0.00217 0.00182 0.00111 0.00192 0.00218 0.00223 ... 
## max values  : 0.019 0.019 0.019 0.020 0.020 0.020 0.019 0.020 0.019 0.019 ...

11.4 Plot

11.4.1 Plot all the area

# Plot Amazon in the year 2020
plot(humidity.rts[['2020']], col=topo.colors(20))

# plot Amazon in October 2020.
plot(humidity.rts[['2020-10-01']], col=topo.colors(20), main="Humidity in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
humidity.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=topo.colors(20), main="Humidity in October 2020")
plot(amaz.basin.shp$geometry, add=T)

11.4.2 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
humidity.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=topo.colors(20))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
humidity.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=topo.colors(20), main="Humidity in October 2020")
plot(amaz.basin.shp$geometry, add=T)

11.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("humidity")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.humidity.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'humidity_working_', '')
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# humidity.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(humidity.freq.na)[3] <- "humidity_na"
# humidity.freq.na <- humidity.freq.na[order(humidity.freq.na$layer)]
humidity.freq.na

11.6 Sample of NA values (February 2002)

11.6.1 Compute

humidity.2002.02.name <- c("Amazon_new_data/7. Specific Humidity/03. Working Data/humidity_working_2002_2.tif")
humidity.2002.02.rast <- rast(humidity.2002.02.name)
humidity.2002.02.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : humidity_working_2002_2.tif 
## name        : humidity_2002_2_proj 
## min value   :          0.004104336 
## max value   :          0.019729620
#
humidity.2002.02.nonNA <- not.na(humidity.2002.02.rast)
humidity.2002.02.nonNA.mask <- mask(humidity.2002.02.nonNA, amaz.basin.shp)
#
humidity.2002.02.nonNA.df <- as.data.frame(humidity.2002.02.nonNA.mask, xy=T)
humidity.2002.02.NA.df <- humidity.2002.02.nonNA.df[humidity.2002.02.nonNA.df[,3]==F,]
humidity.2002.02.NA.df
dim(humidity.2002.02.NA.df)[1]
## [1] 18167

11.6.2 Plot

# Plot
plot(humidity.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))

#
x_min=-0.45e+06; x_max=-0.3e+06; y_min=ymin(humidity.2002.02.nonNA.mask); y_max=1.71e+06
humidity.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#
x_min=0.e+06; x_max=xmax(humidity.2002.02.nonNA.mask); y_min=3.e+06; y_max=ymax(humidity.2002.02.nonNA.mask)
humidity.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#=====================

12 Evapotranspiration data

12.1 Import data

# list of files
amaz.evapotranspiration.list <- list.files("Amazon_new_data/8. Evapotranspiration/03. Working Data",
                                           full.names=TRUE,
                                           pattern = ".tif$")
# Import data with "Terra"
evapotranspiration.rast <- rast(amaz.evapotranspiration.list)
evapotranspiration.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : evapotranspiration_working_2001_1.tif  
##               evapotranspiration_working_2001_10.tif  
##               evapotranspiration_working_2001_11.tif  
##               ... and 237 more source(s)
## names       :  evapo~_proj,   evapo~_proj,   evapo~_proj,   evapo~_proj,   evapo~_proj,   evapo~_proj, ... 
## min values  : 0.000000e+00, -3.037881e-09, -1.953269e-08, -1.696159e-08, -3.509835e-08, -3.830653e-08, ... 
## max values  : 8.157189e-05,  8.477412e-05,  9.205872e-05,  7.918844e-05,  9.691325e-05,  8.340040e-05, ...

12.2 Verification of the data

12.2.1 Verification the Names of layers

names(evapotranspiration.rast) %>% head()
## [1] "evapotranspiration_2001_1_proj"  "evapotranspiration_2001_10_proj"
## [3] "evapotranspiration_2001_11_proj" "evapotranspiration_2001_12_proj"
## [5] "evapotranspiration_2001_2_proj"  "evapotranspiration_2001_3_proj"

12.2.2 Verification of the values

# Verification of the values
evapotranspiration.minmax <- minmax(evapotranspiration.rast) %>% t() %>% as.data.frame()
evapotranspiration.minmax

12.3 Modification of the data

12.3.1 Rename layers

evapotranspiration.rast <- renameLayers(evapotranspiration.rast, 'evapotranspiration_working_', 'evapotranspiration_')
# Order layers
ordered.names.8 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("evapotranspiration_", .)
evapotranspiration.rast <- evapotranspiration.rast[[ordered.names.8]]
# Print
evapotranspiration.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : evapotranspiration_working_2001_1.tif  
##               evapotranspiration_working_2001_2.tif  
##               evapotranspiration_working_2001_3.tif  
##               ... and 237 more source(s)
## names       :  evapo~01_01,   evapo~01_02,   evapo~01_03,   evapo~01_04,   evapo~01_05,   evapo~01_06, ... 
## min values  : 0.000000e+00, -3.509835e-08, -3.830653e-08, -2.597946e-09, -1.043295e-08, -1.211072e-08, ... 
## max values  : 8.157189e-05,  9.691325e-05,  8.340040e-05,  7.563404e-05,  7.022559e-05,  8.013412e-05, ...
cat('\n---------- \n\n')
## 
## ----------
names(evapotranspiration.rast) %>% head()
## [1] "evapotranspiration_2001_01" "evapotranspiration_2001_02"
## [3] "evapotranspiration_2001_03" "evapotranspiration_2001_04"
## [5] "evapotranspiration_2001_05" "evapotranspiration_2001_06"

12.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
evapotranspiration.rts <- rts(evapotranspiration.rast, seq.dates)
evapotranspiration.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  :  0.0e+00 -3.5e-08 -3.8e-08 -2.6e-09 -1.0e-08 -1.2e-08 -5.4e-08  0.0e+00 -1.0e-08 -3.0e-09 ... 
## max values  : 8.2e-05 9.7e-05 8.3e-05 7.6e-05 7.0e-05 8.0e-05 8.3e-05 8.8e-05 7.9e-05 8.5e-05 ...

12.4 Plot

12.4.1 Plot all the area

nbr_colors <- colorRampPalette(brewer.pal(9, "YlOrBr"))
# Plot Amazon in the year 2020
plot(evapotranspiration.rts[['2020']], col=nbr_colors(100))

# plot Amazon in October 2020.
plot(evapotranspiration.rts[['2020-10-01']], col=nbr_colors(100), main="Evapotranspiration in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
evapotranspiration.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=nbr_colors(100), main="Evapotranspiration in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

12.4.2 Plot south of Amazon

nbr_colors <- colorRampPalette(brewer.pal(9, "YlOrBr"))
x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
evapotranspiration.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=nbr_colors(100))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
evapotranspiration.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=nbr_colors(100), main="Evapotranspiration in October 2020")
plot(amaz.basin.shp$geometry, add=T)

12.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("evapotranspiration")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.evapotranspiration.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'evapotranspiration_working_', '')
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# evapotranspiration.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(evapotranspiration.freq.na)[3] <- "evapotranspiration_na"
# evapotranspiration.freq.na <- evapotranspiration.freq.na[order(evapotranspiration.freq.na$layer)]
evapotranspiration.freq.na

12.6 Sample of NA values (February 2002)

12.6.1 Compute

evapotranspi.2002.02.name <- c("Amazon_new_data/8. Evapotranspiration/03. Working Data/evapotranspiration_working_2002_2.tif")
evapotranspi.2002.02.rast <- rast(evapotranspi.2002.02.name)
evapotranspi.2002.02.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : evapotranspiration_working_2002_2.tif 
## name        : evapotranspiration_2002_2_proj 
## min value   :                  -1.071572e-08 
## max value   :                   7.709819e-05
#
evapotranspi.2002.02.nonNA <- not.na(evapotranspi.2002.02.rast)
## 
|---------|---------|---------|---------|
=========================================
                                          
evapotranspi.2002.02.nonNA.mask <- mask(evapotranspi.2002.02.nonNA, amaz.basin.shp)
## 
|---------|---------|---------|---------|
=========================================
                                          
#
evapotranspi.2002.02.nonNA.df <- as.data.frame(evapotranspi.2002.02.nonNA.mask, xy=T)
evapotranspi.2002.02.NA.df <- evapotranspi.2002.02.nonNA.df[evapotranspi.2002.02.nonNA.df[,3]==F,]
evapotranspi.2002.02.NA.df
dim(evapotranspi.2002.02.NA.df)[1]
## [1] 132998

12.6.2 Plot

# Plot
plot(evapotranspi.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))

#=====================

13 Wind Speed data

13.1 Import data

# list of files
amaz.wind.list <- list.files("Amazon_new_data/9. Wind Speed/03. Working Data",
                             full.names=TRUE,
                             pattern = ".tif$")
# Import data with "Terra"
wind.rast <- rast(amaz.wind.list)
wind.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : wind_working_2001_1.tif  
##               wind_working_2001_10.tif  
##               wind_working_2001_11.tif  
##               ... and 237 more source(s)
## names       : wind_~_proj, wind_~_proj, wind_~_proj, wind_~_proj, wind_~_proj, wind_~_proj, ... 
## min values  :    1.069407,    1.048391,    1.137132,    1.105208,    1.024364,    1.063414, ... 
## max values  :    8.202874,    7.747525,    7.947669,    8.901711,    9.390056,    8.901502, ...

13.2 Verification of the data

13.2.1 Verification the Names of layers

names(wind.rast) %>% head()
## [1] "wind_2001_1_proj"  "wind_2001_10_proj" "wind_2001_11_proj"
## [4] "wind_2001_12_proj" "wind_2001_2_proj"  "wind_2001_3_proj"

13.2.2 Verification of the values

# Verification of the values
wind.minmax <- minmax(wind.rast) %>% t() %>% as.data.frame()
wind.minmax

13.3 Modification of the data

13.3.1 Rename layers

wind.rast <- renameLayers(wind.rast, 'wind_working_', 'wind_')
# Order layers
ordered.names.9 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("wind_", .)
wind.rast <- wind.rast[[ordered.names.9]]
# Print
wind.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : wind_working_2001_1.tif  
##               wind_working_2001_2.tif  
##               wind_working_2001_3.tif  
##               ... and 237 more source(s)
## names       : wind_2001_01, wind_2001_02, wind_2001_03, wind_2001_04, wind_2001_05, wind_2001_06, ... 
## min values  :     1.069407,     1.024364,     1.063414,     1.128978,     1.000217,    0.9819801, ... 
## max values  :     8.202874,     9.390056,     8.901502,     8.614891,     7.900918,    7.8636937, ...
cat('\n---------- \n\n')
## 
## ----------
names(wind.rast) %>% head()
## [1] "wind_2001_01" "wind_2001_02" "wind_2001_03" "wind_2001_04" "wind_2001_05"
## [6] "wind_2001_06"

13.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
wind.rts <- rts(wind.rast, seq.dates)
wind.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : 1.07 1.02 1.06 1.13 1.00 0.98 1.05 1.10 1.12 1.05 ... 
## max values  : 8.2 9.4 8.9 8.6 7.9 7.9 8.0 9.1 7.9 7.7 ...

13.4 Plot

13.4.1 Plot all the area

# Plot Amazon in the year 2020
plot(wind.rts[['2020']], col=cm.colors(50))

# plot Amazon in October 2020.
plot(wind.rts[['2020-10-01']], col=cm.colors(50), main="Wind in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
wind.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=cm.colors(50), main="Wind in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

13.4.2 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
wind.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=cm.colors(50))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
wind.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=cm.colors(50), main="Wind in October 2020")
plot(amaz.basin.shp$geometry, add=T)

13.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("wind")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.wind.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'wind_working_', '')
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# wind.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(wind.freq.na)[3] <- "wind_na"
# wind.freq.na <- wind.freq.na[order(wind.freq.na$layer)]
wind.freq.na

13.6 Sample of NA values (February 2002)

13.6.1 Compute

wind.2002.02.name <- c("Amazon_new_data/9. Wind Speed/03. Working Data/wind_working_2002_2.tif")
wind.2002.02.rast <- rast(wind.2002.02.name)
wind.2002.02.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : wind_working_2002_2.tif 
## name        : wind_2002_2_proj 
## min value   :         1.199288 
## max value   :         9.328015
#
wind.2002.02.nonNA <- not.na(wind.2002.02.rast)
## 
|---------|---------|---------|---------|
=========================================
                                          
wind.2002.02.nonNA.mask <- mask(wind.2002.02.nonNA, amaz.basin.shp)
## 
|---------|---------|---------|---------|
=========================================
                                          
#
wind.2002.02.nonNA.df <- as.data.frame(wind.2002.02.nonNA.mask, xy=T)
wind.2002.02.NA.df <- wind.2002.02.nonNA.df[wind.2002.02.nonNA.df[,3]==F,]
wind.2002.02.NA.df
dim(wind.2002.02.NA.df)[1]
## [1] 18167

13.6.2 Plot

# Plot
plot(wind.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))

#
x_min=0.e+06; x_max=xmax(wind.2002.02.nonNA.mask); y_min=3.e+06; y_max=ymax(wind.2002.02.nonNA.mask)
wind.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#=====================

14 Air Temperature data

14.1 Import data

# list of files
amaz.airtemp.list <- list.files("Amazon_new_data/10. Air Temperature/03. Working Data",
                                full.names=TRUE,
                                pattern = ".tif$")
# Import data with "Terra"
airtemp.rast <- rast(amaz.airtemp.list)
airtemp.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : airtemp_working_2001_1.tif  
##               airtemp_working_2001_10.tif  
##               airtemp_working_2001_11.tif  
##               ... and 237 more source(s)
## names       : airte~_proj, airte~_proj, airte~_proj, airte~_proj, airte~_proj, airte~_proj, ... 
## min values  :    268.7742,    270.8836,    270.5267,    270.5608,    269.2101,    269.1326, ... 
## max values  :    303.1154,    304.1908,    304.2381,    303.4955,    304.7826,    304.7043, ...

14.2 Verification of the data

14.2.1 Verification the Names of layers

names(airtemp.rast) %>% head()
## [1] "airtemp_2001_1_proj"  "airtemp_2001_10_proj" "airtemp_2001_11_proj"
## [4] "airtemp_2001_12_proj" "airtemp_2001_2_proj"  "airtemp_2001_3_proj"

14.2.2 Verification of the values

# Verification of the values
airtemp.minmax <- minmax(airtemp.rast) %>% t() %>% as.data.frame()
airtemp.minmax

14.3 Modification of the data

14.3.1 Rename layers

airtemp.rast <- renameLayers(airtemp.rast, 'airtemp_working_', 'airtemp_')
# Order layers
ordered.names.10 <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% 
  format(., '%Y_%m') %>% 
  paste0("airtemp_", .)
airtemp.rast <- airtemp.rast[[ordered.names.10]]
# Print
airtemp.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## sources     : airtemp_working_2001_1.tif  
##               airtemp_working_2001_2.tif  
##               airtemp_working_2001_3.tif  
##               ... and 237 more source(s)
## names       : airte~01_01, airte~01_02, airte~01_03, airte~01_04, airte~01_05, airte~01_06, ... 
## min values  :    268.7742,    269.2101,    269.1326,    270.1584,    269.9435,    269.3073, ... 
## max values  :    303.1154,    304.7826,    304.7043,    304.4511,    303.5922,    304.0975, ...
cat('\n---------- \n\n')
## 
## ----------
names(airtemp.rast) %>% head()
## [1] "airtemp_2001_01" "airtemp_2001_02" "airtemp_2001_03" "airtemp_2001_04"
## [5] "airtemp_2001_05" "airtemp_2001_06"

14.3.2 Create ‘rts’ object

# Create a sequence date for 'rts' object
seq.dates <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month")
airtemp.rts <- rts(airtemp.rast, seq.dates)
airtemp.rts
## Raster Time Series with monthly periodicity from 2001-01-01 to 2020-12-01 
## class       : SpatRasterTS 
## raster filename    : /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2001_1.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2001_2.tif, /Users/ABIDM/Documents/Project/Code_official_data/Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2001_3.tif ... 
## raster dimensions  : 5860, 7806, 240  (nrow, ncol, nlyr)
## raster resolution  : 500, 500  (x, y)
## raster extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## min values  : 269 269 269 270 270 269 269 269 270 271 ... 
## max values  : 303 305 305 304 304 304 304 305 305 304 ...

14.4 Plot

14.4.1 Plot all the area

# Plot Amazon in the year 2020
plot(airtemp.rts[['2020']], col=colorRamps::matlab.like(100))

# plot Amazon in October 2020.
plot(airtemp.rts[['2020-10-01']], col=colorRamps::matlab.like(100), main="Air Temp. in October 2020")
plot(amaz.basin.shp$geometry, add=T)

#
airtemp.rts[['2020-10-01']] %>% 
  mask(mask = amaz.basin.shp) %>% 
  plot(col=colorRamps::matlab.like(100), main="Air Temp. in October 2020")
## 
|---------|---------|---------|---------|
=========================================
                                          
plot(amaz.basin.shp$geometry, add=T)

14.4.2 Plot south of Amazon

x_min=-1.25e+06; x_max=0.05e+06; y_min=ymin(burntArea.rast); y_max=2.4e+06

# Plot south of Amazon in the year 2020.
airtemp.rts[['2020']]@raster %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::matlab.like(100))
## 
|---------|---------|---------|---------|
=========================================
                                          

|---------|---------|---------|---------|
=========================================
                                          

# Plot south of Amazon in October 2020.
airtemp.rts[['2020-10-01']] %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  mask(mask = amaz.basin.shp) %>%
  plot(col=colorRamps::matlab.like(100), main="Air Temp. in October 2020")
plot(amaz.basin.shp$geometry, add=T)

14.5 Verification of the NA values

# cl <- makeCluster(detectCores() - 1)
# registerDoParallel(cl, cores=detectCores() - 1)
# 
# tic("airtemp")
# #Raster to datatable in parallel: one raster per thread
# rasList <- foreach (ras_id=amaz.airtemp.list, .packages=c('terra', 'sf'), .combine='c') %dopar% {
#   #Read all rasters into one big stack
#   ras <- rast(ras_id)
#   # Rename layers
#   ras <- renameLayers(ras, 'airtemp_working_', '')
#   #
#   ras.nonNA <- not.na(ras) 
#   ras.nonNA.mask <- mask(ras.nonNA, amaz.basin.shp)
#   ras.freq.na <- freq(ras.nonNA.mask, digits=0, value=0, usenames=T)
#   
#   list(ras.freq.na)
# }
# stopCluster(cl)
# toc()
# 
# #Bind all per-raster into one dataframe
# airtemp.freq.na <- rbindlist(rasList, fill=T, use.names=T)
# #
# colnames(airtemp.freq.na)[3] <- "airtemp_na"
# airtemp.freq.na <- airtemp.freq.na[order(airtemp.freq.na$layer)]
airtemp.freq.na

14.6 Sample of NA values (February 2002)

14.6.1 Compute

airtemp.2002.02.name <- c("Amazon_new_data/10. Air Temperature/03. Working Data/airtemp_working_2002_2.tif")
airtemp.2002.02.rast <- rast(airtemp.2002.02.name)
airtemp.2002.02.rast
## class       : SpatRaster 
## dimensions  : 5860, 7806, 1  (nrow, ncol, nlyr)
## resolution  : 500, 500  (x, y)
## extent      : -2156811, 1746189, 1625314, 4555314  (xmin, xmax, ymin, ymax)
## coord. ref. : South_America_Albers_Equal_Area_Conic 
## source      : airtemp_working_2002_2.tif 
## name        : airtemp_2002_2_proj 
## min value   :            269.4591 
## max value   :            304.5671
#
airtemp.2002.02.nonNA <- not.na(airtemp.2002.02.rast)
## 
|---------|---------|---------|---------|
=========================================
                                          
airtemp.2002.02.nonNA.mask <- mask(airtemp.2002.02.nonNA, amaz.basin.shp)
## 
|---------|---------|---------|---------|
=========================================
                                          
#
airtemp.2002.02.nonNA.df <- as.data.frame(airtemp.2002.02.nonNA.mask, xy=T)
airtemp.2002.02.NA.df <- airtemp.2002.02.nonNA.df[airtemp.2002.02.nonNA.df[,3]==F,]
airtemp.2002.02.NA.df
dim(airtemp.2002.02.NA.df)[1]
## [1] 18167

14.6.2 Plot

# Plot
plot(airtemp.2002.02.nonNA.mask, col=c("darkorange1", "forestgreen"))

#
x_min=0.e+06; x_max=xmax(airtemp.2002.02.nonNA.mask); y_min=3.e+06; y_max=ymax(airtemp.2002.02.nonNA.mask)
airtemp.2002.02.nonNA.mask %>% 
  crop(ext(x_min, x_max, y_min, y_max)) %>% 
  plot(col=c("darkorange1", "forestgreen"))

#=====================

15 Merge all dataframes of NA values.

# create the dataframe
seq.dates.str <- seq(as.Date("2001-1-1"), as.Date("2020-12-1"), by = "month") %>% format("%Y_%m")
amaz.na.df <- as.data.frame(seq.dates.str)
colnames(amaz.na.df) <- "layer"
# Merge
amaz.na.df <- list(amaz.na.df, 
                   burntArea.freq.na[,-2], 
                   landCover.freq.na[,-2], 
                   precipitation.freq.na[,-2],
                   soilmoisture.freq.na[,-2],
                   elevation.freq.na[,-2],
                   landsurftemp.freq.na[,-2], 
                   humidity.freq.na[,-2], 
                   evapotranspiration.freq.na[,-2],
                   wind.freq.na[,-2], 
                   airtemp.freq.na[,-2]) %>% 
  reduce(full_join, by="layer")

amaz.na.df
amaz.na.df[133:144,]

#=====================